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ABSTRACT 

Three statistical theories of turbulence are applied to thermal convection 
between infinite slippery plates in the limit of infinite Prandtl number. The 
range of Rayleigh number R investigated is 657 <R^1.25x 10 4 . The theories, 
which are compared in a series of numerical calculations, are the direct inter- 
action approximation, the quasi-linear approximation, and the quasi-normal ap- 
proximation. The direct interaction approximation gives results for the evolu- 
tion Nusselt number in good agreement with some simple exact solutions to the 
problem. The flow predicted by this method is very presistent. Turbulence 
first appears just at critical R and its intensity gradually increases with increas- 
ing R. The quasi-normal approximation gives satisfactory results for R <2000, 
but some of the initial value problems lead to an unphysical negative tempera- 
ture autocorrelation spectrum for larger R. In none of the initial value problems 
for R > 2 x 10 3 did the quasi-normal procedure yield a sensible stationary state. 
The quasi-linear approximation appears to upperbound the heat transport given 
by both the other approximations predicting about 10% larger heat transport than 
the direct interaction approximations. 
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I. INTRODUCTION 


This paper presents results of an application of the direct interaction (DI) 
approximation'*’ to thermal convection at infinite Prandtl number, in a fluid con- 
fined by infinite, slippery, infinitely conductive horizontal boundaries. Results 
are presented for initial value problems in which the temperature and velocity 
covariances are allowed to evolve to their steady state values. The infinite 
Prandtl number (o — °°) regime is particularly simple because the non-linear terms 
in the momentum equation vanish, so that one has to study a single non-linear, 
scalar equation for the covariance of the temperature field. Since the direct 
interaction equations are formidable, it is attractive to study first a relatively 
tractable situation to assess the method’s realism before addressing the more 
interesting and difficult case of convection at finite Prandtl number o-. 

There is also presented a comparison of the DI approximation with the 
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"quasi-linear" approximation, the "quasi -normal" approximation, and a direct 
integration of the amplitude equations for a reasonably large ensemble of initial 
data. The first three approaches are methods of closing the hierarchy of moment 
equations. 

The quasi-linear or "mean-field" approximation, discards third order cumu- 

lants. This is equivalent to neglecting deviations of the bilinear terms in theNavier 

Stokes equations from their horizontal averages. The method appears to have 

2 4 

some quantitative validity, for high Prandtl number fluids. * It is known to 
have several defects, including the inability to assess the importance of the 
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terms it omits. For systems which have no mean fields, this method gives the 

trivial answer of dissipative decay of the initial flow and temperature fields. 

3 

The quasi-normal approximation evaluates fourth order moments as if the 
flow and temperature fields were normally distributed. This approach has been 

5 

investigated for isotropic turbulence by Ogura, who solved some of the initial • 

value problems and found the unphysical result of a partly negative kinetic energy 
spectrum. However, failure for the isotropic problem does not preclude success 
in problems in which there is a mean field present. Interaction of fluctuations 
with the mean field may stabilize the system. 

The direct interaction approximation may be viewed as a logical improve- 
ment on the quasi-normal approximation. Both approximations give quantitative 
values for an effective dissipation acting on each degree of freedom of the fluid. 

The expressions for this dissipation involve functions which represent the average 
response of a mode to an infinitesimal initial perturbation. The quasi-normal 
approximation yields response functions determined wholly by molecular dissi- 
pation. The direct interaction approximates the effects of eddy viscosity, in 
a self consistent fashion. 

The direct numerical integrations of the equations of motion serve to assess 
the validity of the statistical methods described in the preceding paragraphs, # 

especially in lieu of any experimental data for slip boundary conditions. These 

» 

results were obtained by constructing a random set of initial data, allowing each 
realization of the flow to evolve in three dimensions, and then ensemble— averaging 
the results, thereby obtaining values for the heat transport, and correlation functions. 
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The numerical results suggest that the direct interaction approximation 


may be satisfactory for thermal turbulence. The convective heat transport 
predicted is in good agreement with the direct integrations of the Boussinesq 
equations. The auto-correlation functions have very large relaxation times, 
suggesting an almost static flow for the range of Rayleigh numbers investigated 
(R< 1.25 x 10 4 ). The behavior of the steady-state Green's functions indicate that 
the system is close to being marginally stable. That is to say, a large-scale 
temperature perturbation decays very slowly under the joint action of the mean 
field and the eddy-dissipation processes. The fact makes the predictions of the 
direct interaction approximation agree closely with those of the quasilinear 
method with regard to the steady state value of the Nusselt number. The Green's 
functions predicted by the two methods are completely different. 

The quasi-normal approximation appears to be a reasonable approximation 
only for R < 2000. For R > 2000, the approximation predicts the development of 
a negative value for the square of the temperature fluctuation field for some 
initial data investigated. It thus appears that the stabilizing influence of the 
fluctuation-field-mean-field interaction is very limited, and the eddy-type terms 
introduced by the direct interaction represent an essential feature. We should 
point out, however, that our results for the initial value problem do not preclude 
the validity of ‘the quasi-normal approach for steady-state turbulence at any 
value of R. 
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H. THE EQUATIONS OF MOTION 


Under the stated boundary conditions, the infinite -cr Boussinesq equations 
for the horizontally averaged temperature field T(z , t ) and deviation 6 ( r , t ) of 
temperature T(r\ t) from T( z , t) may be reduced to 



n, a 


where; 


J n a - (R/tt 4 ) a 2 /(n 2 + a 2 ) 2 


B 


n m 



n - m j 



/3 - nwT 

n n 


and 


M aa'a" = £ aa"a" g + F “a'a" /g + g 

npq pq n,p + q pq \ p,n+q q, n+p 


with 
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In Eqs. (1) and (2), 6 n a is the component of 8(r, t) proportional to e ia ' x sin nnz t 
and T n is the component of T(z, t) proportional to sin7rz . A brief derivation 
of Eqs. (1) and (2) from the Boussinesq form of the Navier Stokes equations is 
given in Appendix A. 

The nonlinear terms in (1) are sorted into two groups: the terms propor- 
tional to B nm , called the mean field terms, and the terms proportional to M, 
called the fluctuating self interactions. The M coefficients satisfy 

M aa ' a " + M a ' a ~ a " + M a " a_a ' = 0 , 

npq pnq qnp 

from which it follows that the volume -averaged square of the temperature field 
satisfies 


+ + 2> 0|tJ 

n, a n, a n 


1+7 T 


L ■*. 


( 3 ) 


Equation (3) states that the flux of heat into the system (the right hand side) is 
equal to the time rate of accumulation of the heat energy plus the rate of 
molecular dissipation of macroscopic motion into heat. 
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HI. STATISTICAL THEORIES OF TURBULENT CONVECTION 


Our primary goal is to determine the temperature autocorrelation function 
matrix, 


C.('l*') = i e ° <*) e „"“ (•')> ' (4) 

The average ( ) is over an initial ensemble, in which the 9 * are assigned normally 
distributed values, and have a covariance matrix (0 | 0). Statistical homo- 
genity in the horizontal, and identity of horizontal averages to ensemble averages 
are assumed. Homogeneity in the horizontal implies 

(6> n a (t ) 6> m a ' (t')) = 0 , unless a + a' = 0 . 

The heat transport per unit area through the convecting layer is determined by 
(2) and the assumed equivalence between horizontal and ensemble averages. 

An equation for > p may be obtained by multiplying (1) by 9~ a ( t ' ) and averaging: 

(ar + ('I o = T V B „ P I »'> 


+ 


L 

p, q 


+ a' 


M aa ' a " 3 aa ' a " ■ ( t | t ' ) , 

npq npq v 1 7 


( 5 ) 


where 


I f) = (t)). 
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The remainder of this section outlines several theoretical approaches for 


obtaining closed equations for expressing the triple moment 0 as a functional of 

i p and B. The particular closure approximations to be studied are the quasi- 

linear approximation (in which 0 = 0), quasi-normal approximation, and the 

direct interaction approximation. All three of these approximations have been 
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discussed in some detail elsewhere, * * and the discussion to follow is in- 
tended to form a unified reference frame for the remainder of the paper. The 
only new possible ingredient is in the discussion of the quasi-normal approxima- 
tion. Here, some new non-diagonal correlation functions \fj nm occur, which are 
zero in isotropic turbulence. 

a. The Quasi-Linear Approximation 

If the horizontal exchange of heat and energy is very weak compared to the 

vertical exchange induced by the term B0 it may be argued that 3 - 0 is a good 

approximation. This procedure — called hereafter the quasi-linear approximation — 

2 4 

has been investigated in some detail. * It appears to be a fairly accurate and 
consistent method for large Prandtl number fluids. The mean-field terms, in 
an initial value problem, inhibit the continued growth of an initial perturbation 
by modifying the mean field through (2). This feedback mechanism stabilizes 
the system, yielding an eventual steady state. It is not clear that once the 3 term 
is reintroduced, the mean-field terms still play this same role. 
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b. The Quasi-Normal Approximation 


A systematic formal procedure for including approximations to the triple- 

0 

moment transfer term is discard of cumulants procedure. An hierarchy of 
equations for all the moments with simultaneous time arguments is first obtained. 
An equation for 3 is obtained by multiplying~(l), by 6~ a ' ( t ) 6~ a " ( t ) , en- 
semble averaging, and adding the result to similar expressions obtained from 
(1) by permuting the indices n, p, and q. The equation for 3 thereby obtained is 



+ v a 

p 



aa ' a " 

° mpq 


D g a rj aa ' a" | ga'g aa'a" 
mr'rpq pr C mrq 


+ B a "3 aa ' a 

qr mp r 


) 



/ V / III t V II III I V 

+ M a a a 0~ aa a a 

prv Xnqrv 


M a a 

qr v 


^ -aa a a 
J mprv 




where 


q aa ' a " a 1 " = /q a Q a ' Q a " Q a"'\ 

y nmp r \ n m p r / 

is a fourth order moment. An equation for it may be obtained by an algorithm 
like that used to obtain the 3 equation. 

By proceeding in this way, a hierarchy of equations is obtained which re- 
late the moments of order n to those of order n + 1. To obtain closed equations 
of motion for a given moment (for example 3) the higher-order moment (Q) is 
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eliminated from the equation of motion by assuming that the form of the higher 
order moment is the same as it would be if the <9 n a (t) were Gaussianly dis- 
tributed. Thus, in the fourth order closure. 


/ n in 
q aa a a 

^nmpq 


dj a \b a " 8 , 8 « ,„ + 0 a l // a ' S „ 8 , 

r nm r a, -a a ,~a ’‘np’mr a, -a a ,-a 


+ a ' <Jj a ‘ 8 in 8 i i, . /rr\ 
f nr y mp a, “a a ,-a (I) 


Equations (5), (6), and (7) now comprise a complete set of equations for the de- 
termination of \p( t , t ) . 

The closure procedure using (6) and (7) is a generalization of the quasi- 
normal approximation — originally proposed for isotropic turbulence by Tatsumi 
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and by Proudman and Reid. Ogura has shown numerically that the method is 
not satisfactory for isotropic turbulence; an integration of some sample initial 
data developed an energy spectrum with negative regions after a short time, for 
moderate Reynolds numbers. 

Higher-order cumulant discard approximations may in principle be an im- 
provement over the Tatsumi-Proudman-Reid quasi-normal approximation, but 
to integrate numerically any such higher order scheme appears prohibitively 
difficult. 

The quasi-normal equation for 0 may be put into more convenient form by 

O 

introducing the Green's function G nm ( t, t ' ), which satisfies the following 
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equations: 


(£ + -„ a ) G„“. 


(t I t') 


L 


B“G ( t t ' ) , 

np pm v 1 ' 


with 


G (t t) 

n m v ’ / 


(8) 


Using (6), (7), and (8), 3 may be evaluated as 


rr a a a 

C npq 




( 1 1 1 ' ) g ( 1 1 1 ' ) G q a ; ( 1 1 1 * ) s a a b a ; a " ( t - > 


where; 


c aa a 
ab c 


T fM; r a „ a ' a " + ^, a ' 0 a " 

/ \arv avv / ^br ^cv 


+ 2 7 /M, a '“ aa " + M L a ' a " _a ^ 0 a 0 a " . (9) 

/ \ b r v bvr y “a r ”c v 

c. The Direct Interaction Approximation 

In the DI approximation, the evaluation of 3 is accomplished through the 
intermediary of a Green's function, g n a m ( 1 1 t ' ), whose role in the theory is 

o 

similar to the G function introduced in the quasi-normal theory. This function is 
defined as the ensemble-averaged response of mode 0 * at time t to an infinitesimal 
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impulsive perturbation in mode 0 m a at time t ' . In terms of the g*s, the direct 

rr 7 

interaction approximation for 0 is 


3 | f) 


ft' 

M; r V a '" a " J dt« g m a s ( t ' | t")^ r ' (t | t")^ q a v "(t I t") 

s, r,v ® 


L "•"*■*' f 

s , r , v 


dt"g“ (t | | (t t") 


L f 

s, r , v 


The DI approximation for g n a m (t | t ' ) is 


dt" gq * v (t | t")0 m * r (t' | t")0 p ‘ (t | t") ( iQa) 


^ ,_8 





+ 



s, r, v 


a = a ' + a " 


t')g q a s "(t I t")0 p v (t 1 1") 


(10b) 


with g n a m (t I t) - 5 m , 

Equations (10a) and (10b) determine i p* m (t | t') given the initial values of 
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The DI equation for g n ^ bear some similarity to that for G nm . Equation (10b) 
contains an additional convolutional term which represents the relaxational effect 
of eddy dissipation. In the limit R -0, g n a m and G n “ become identical. 

The structures of the quasi-normal and direct interaction approximations 
are similar with regard to 3(t | t') also. The quasi-normal equations may be 
obtained from the DI equations by making on the latter the following alterations: 

1. Change g n “ ( 1 1 t ' ) to G n “ ( t | t ' ) by deleting the last term of (10b) . Then, 

2. Replace 0 n “ ( t | t ' ) in (10a) according to 

Cm ^ I *'> ^ ^ ^ 1 t ' ) ^ (t ' ) • (11) 

a 

In the direct interaction equations, (5), (10a), and (10b), three basic terms 
determine the behavior of 0 and g: (1) the dissipation terms, proportional to 
v n’ (2) the mean-field terms — the first terms on the right hand side of all the 
equations, and (3) the fluctuating self-interactions, which are the last terms in 
the equations. Only the first two of these terms enter the net entropy balance. 

A physical interpretation of the direct-interaction treatment of the ^-self- 
interactions may be given in terms of a generalization of eddy-dissipation con- 
cepts. The right hand side of (10a) consists of basically two types of terms: 

(1) an input term to the mode a, and two drain terms from a. The input terms 
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contain the factor 


■ •.* ' 



t " ) dt 


and the drain terms contain the factors 



t " ) dt * 


The conservative nature of all three statistical procedures discussed here 
may be verified directly from (5), (9) or (10a), on using the identity the M coef- 
ficients satisfy. For the direct-interaction method the conservation law is in- 
dependent of the values of g n a m . 

IV. SOME GENERAL COMMENTS ON THE. STATISTICAL THEORIES 

Before proceeding to the numerical results, it is worthwhile to consider 
first certain consistency properties and some general properties of the solu- 
tions for the statistical approximations. 

The most obvious consistency property to examine is the conservation of 
entropy (3). All the approximation considered here satisfy this equation, as 
stated previously. A second, and perhaps less obvious consistency condition on 
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the methods has to do with bounds on the >/' nm ( t | t ' ) , which arises from the 
fact that these quantities are ensemble averages of products of real temperature 
fluctuation fields, 6 ( r , t ). From this fact and by Schwartz 1 inequality it follows 
that 


1 Vtn (^’1 t ) ^ ( t 'Tt 7 ) 


or defining the correlation coefficient. 


R 


n m 


(t I t') 


(t 1 t ') 

^ „„(tl ft T > 


KnmCH t') < 1- 


For the quasi-linear approximation, R nm ( t ) t ' ) — 1 because 0 is an exact 
solution of amplitude equations, without phase mixing (fluctuating self interactions 
are discarded). 

The quasi-normal approximation deals only with the time -diagonal terms 
Aim (tit)- It does not assure R nm (t | t) < 1. 

In the direct interaction approximation, R nm ( t j t ' ) < 1 for all t , t'. This 
follows from the fact that the approximation corresponds to an exact solution for 
a set of model amplitude equations. This point has been extensively discussed 

9 

by Kraichnan. 

The coefficients R iim ( t | t') provide a measure of the rate of change of the 
temperature and velocity patterns in a typical realization. Consider a statistically 
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stationary situation, for which R nm ( t | t ' ) depends only on the time difference 
r - 1 1 - t ' |. If R nm (r) decays rapidly as r increases, this rate of change is 
large. On the other hand, if R nra (r) decays very slowly, there is nearly steady 
cellular convection. If R nm is independent of r, the system is said to be static. 

A static form for R hm need not preclude a random irregularity of the spatial 
pattern of the flow. 

The direct interaction approximation may admit of static solutions, as has 
been pointed out by Kraichnan. 1 ^ To see that the assumption of t- and t '-independent 
0 (t 1 1 ' ) is consistent with the direct interaction equations in the limit t -oo f 
it is only necessary to examine (5), (10a), and (10b) and recall that as t-oo both 
4> and g should relax to functions of argument r = t - t' only. Then (5) and (10a) 
imply that the static solution for \p is a function of the integral of g, 

§nm = f ds gnm( S ) ds ' 

Jo 

An equation for § nra may be obtained by integrating (10b) over t, and using the 
boundary conditions on g. This equation is 



P 


+ 


L 


M nps M s 


§ a 

vm 


e a 0 a 

^qs r pr 


( 12 ) 
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assuming, of course that the integral of g exists. The static solution is thereby 
reduced to a set of algebraic equations. The question of whether these solutions 
are stable as compared to the dynamic solutions is unanswered. Conceivably the 
static solutions are stable below a certain critical Rayleigh number while the 
dynamic solutions are stable above. 

The quasi-normal approximation cannot distinguish static and dynamic 
statistically stationary flows because it deals with only simultaneous correlation 
coefficients R nm (t | t). 

A third consistency requirement comes from the requirement that the tem- 
perature field T(r, t) be limited by the boundary temperatures T(0, t) = 0, and 
T( 1, t ) = -1.0, provided the initial temperature field is so limited. This comes 
about because there are no heat sources in the fluid. This constraint on T( r, t) 
implies 


0 < <([" T(7, t)] n ) < l(n = 1, 2, •••) 

If the constraints are not imposed at t = 0, they are applicable only as t -oo > 
in which limit the smoothing action of the diffusivity erases transient behavior 
making the statistical moments settle down to final state values. 

None of the statistical approximations discussed in this paper give a priori 
assurance of conforming to all the moment bounds onT. The only method for 
which results are available is the quasi-linear approximation, and they indicate 
violations of the bounds for certain problems. Indeed, the contrary would be 
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surprising, since this approximation omits some of the terms responsible for 

the bounds. Recently, Durney^ 1 has applied the quasi-linear approximation to a 

study of convection in spherical shells. He finds that for R£ 2.4 x 10 4 the odd- 

ordered moment bounds are violated in a thin region in near the boundary layer. 

The violation is slight but still none the less present. 

For the problem of plane parallel convection the quasi-linear approximation 

does not lead to violation of the bound on T for either free or rigid boundaries, 

or for any value of R, provided the system has reached the steady state. This 

12 

statement is in apparent contradiction to conclusions reached by Deardorff, 
who reported violations of the bound conditions near the boundary, for a par- 
ticular initial value integration. 

V. RESULTS 

a. Numerical Procedure 

The numerical procedure for integrating the statistical equations for 
</> (t 1 1') is described in detail in Appendix C. The complexity of the equa- 
tions imposes severe restrictions on the number of modes explicitly treated, 
and on the range of Rayleigh number, R. The present paper restricts itself to 
only five vertical wave numbers, and up to three horizontal wave numbers to 
describe the \b a ft I t' Afield. The results of the next two sub-sections are 
obtained by approximating the a-dependence of «//(a) by a constant value for 
1//2 < a < V2 and zero elsewhere. In section Vc, in which greater accuracy is 
desired three values of a are included to represent </'(a). Details of both methods 
of numerical approximation are given in appendix C. 
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to difference the equations, and the trapizoidal rule to evaluate the convolution 
integrals in (10a), and (10b). The procedures guarantee exact conservation 
properties for the system. 

b. Results for the Direct Interaction Method 

The first point to establish concerning the validity of the direct interaction 
method in whether or not it faithfully reproduces the onset of convection just 
above the critical Rayleigh number, R c = 657.71. To test this point, three runs 
were made at R = 700, 830, 1000, using only two vertical modes and a single a. 
These runs used unaveraged J n a , and v ^ coefficients with a = 1//2, in order to 
facilitate a comparison of the numerical results with results known from 
perturbation theory. The parameters of these runs were a = 1//2, a l = V2, 
and A = 0.20. The initial value the </» nm spectrum was 

= ( 1_R C / R ) (13) 

0 = 0 , for n / m . 

The initial T n field was such that (dT n /dt) = 0 for a ^ nm field given by (13). The 

above </> spectrum is that implied by perturbation theory, as for example that of 

13 

Malkus and Veronis. The values of N u for t = 8.0 are shown in Fig. 1, and are 
compared with results of other statistical theories. This value of t is sufficient 
in these cases to assure l'/ , n /</' n l - 0.006. 

The values of R u (t | t') and R 22 (t | t') for the R = 10 3 run is shown in 
Fig 2. These curves indicate very long correlation times, which implies a very 
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presistent flow. The Green's functions g n (t | t ') and g 22 (t | t' ) are shown 
in Fig. 3. These results are in good agreement with the analysis of the two mode 
static system analyzed in Appendix B (cf. B12 and B13), which indicate an 
extremely slow decay of small perturbations introduced in the flow. Results for 
the R(t | t ') and g(t | t') at the other values of R show an entirely similar 
behavior. 

Some results for R = 4 x 10 3 are shown in Figs. 4-7. The parameters, for 
this case are, M = 5, a Q = 1//2, <x i = /2, and A = 0.05. The initial value for 
i/» is that obtained from the symmetrized single - a hexagon solution. Figures 
4, and 5 gives the correlation functions R^ (t I t'), for t = 1.5. We surmise 
from these figures that at R - 4 x 10 3 ,the solutions are no longer nearly static 
in character. This point is not entirely clear from the numerical analysis 
since the solutions have not entirely reached steady state at t = 1.5, and con- 
tinuing the calculation beyond t = 2 by the time forward integration technique 
becomes prohibitive in terms of machine time. 

The Green's functions are given in Figs. 6, and 7. Again, note the slow 
relaxation of the large scale mode, g n . A striking feature seen in these curves 
is the qualitatively different behavior between the even-even modes and the odd- 
odd modes. All the odd-odd Green's functions are "locked in" to the behavior 
of G 1 and have quite long relaxation times. None of them develop a negative 
value. The even-even modes are also "locked'in" to the behavior of G 22 . They 
all develop negative values, and their time integral is quite small compared, say, 
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to V-\ This behavior to reminiscent of that predicted by (B12) and (B13), and 
makes the even-even <, p nm (t | t) correlation matrix elements very small. 

Figures 8 and 9 give some results for R = 10 4 . These results are similar to 
the corresponding ones found at R = 4 x 10 3 , but the system appears to be less 
static. This is suggested by comparing the correlation coefficients as given in 
Figs. 4 and 9. 

c. Results for the Quasi-Normal Method 

The equations for this method are obtained from (5), (10a), and (10b) by 
eliminating i/>(t | t') using (11), and replacing the G equations by (8). Numerical 
results are first presented for the two mode system described at the beginning 
of Section V (b) . 

For R< 2000 the two-mode quasi-normal results appear to be satisfactory. 
The evolution of the system from the initial state prescribed in (12) is reasonable 
and closely resembles the direct interaction results. Values for the Nusselt 
number, N u , are given in Fig. 1, and are compared to results for the other 
methods investigated. The steady state values of N u are consistently lower for 
this method than for the DI or quasi-linear method. 

For r£ 2000 the quasi-normal procedure does not give sensible results. 

The predicted behavior of 0 n (t ) is shown in Fig. 10 for R = 2000. It is seen 
thati/^ (t) evolves into an unphysical -negative region after t- 2.0. The 
numerical accuracy of this curve is not too great especially near the point at 
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which i p n - 0. The time step size used in this graph is A = 0.15. It was found 
that decreasing the step size made the behavior of \p near t = 2.0 more pre- 
cipitous, so that the actual long time behavior of < p n (t ) is slightly more 
singular than is shown in Fig. 10. 

One source of this trouble is not difficult to trace. Recall first that the 
value of N u for this method is below the quasi-linear result, for which v x - B u . 

o 

This means that for the steady state G n (t), as given by (8) has an ex- 
ponential growth rate X = -v l + B n > 0. The numerical results show, 
reasonably enough that X increases withR. Now look at (Bl) as it 
specializes to the quasi-normal case. The net contribution of the fluctuating 
terms to is, 

2M U2 f SiCt I s)gi(t: 

Jo 

As the steady state is approached, the Green's function product becomes 

(2X-v )(t-s) 

e .If now \p 2 > as is only sensible at small R, the above term 

will remain finite as t only if 


s)i 2 (t I s) (<// 1 (s) 0 2 (s) (s)0 x ( S ) ) 


x < ^2/2 • 
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Thus, as R is increased above the value for which A. = v 2/ /2, either \p n (t | t ) 
will eventually become negative as t increases or \p 2 > </>j. The numerical re- 
sults indicate this value of R to be —2,300. 

These simple results suggest that the quasi-normal method suffers yet 

5 

another defect than that pointed out by Ogura; an inability to cope with systems 
having driving forces produced by an instability in the flow field. Of course, it 
may be argued that these results do not survive if a more realistic flow field 
consisting of many Fourier modes is studied. This does not appear to be the 
case, as is discussed in section Vd. 

d . Numerical Integration of the Amplitude Equations 

Although experimental results are not yet available for free boundaries, the 
results of the statistical methods can be compared to the numerical integrations 
of the equations for T(r, t). 

The problem to be considered is: given an initial ensemble of amplitudes 
T n “ i (0) (where i = (1, • • • , N) labels the ensemble member) determine the 
ensemble average covariance, 

N 

Cel' 1 ) 5 w <»>"." <«’> - 

i= l 

Here (> is the fluctuation of the temperature field T from the horizontally averaged 
value, T n The equations of motion for 0 n a . and T n . are (1) and (2). The initial 
data for the complex amplitude iV a . are Gaussianly distributed in the index (i) 
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with a specified covariance matrix i/^ (0). The initial values for T n . is such 
that (dT n .y'dt ) =0. These initial data are consistent with those used for 6 and 
T in the statistical treatment of the problem. 

The correspondance between amplitude and statistical calculations is com- 
plete only if the ensemble is homogenous and isotropic in the horizontal. This 
condition permits a replacement of horizontal averages by ensemble averages, 
and is useful in closing the statistical hierarchy of moment equations. By ap- 
pealing to this equivalence <9 may be interpreted as a fluctuation from an ensemble 
average . 

The 6 n a . were selected from a Gaussianly distributed complex set of numbers 
(with real covariance 0 nm (0)) constructed by first generating a uniform distri- 
bution of numbers, and using this uniform set in connection with the central limit 
theorem to obtain the Gaussian distribution. Care was taken to assure statistical 
independence of real and imaginary parts of 9 . 

The condition of isotropy in the horizontal may be realized in a numerical 
experiment only if the number of a_’s included is very large. Since machine 
calculations necessarily treat a finite number of a_’s it is important to have 
some measure of the degree of isotropy of the experiment. We shall measure 
this by a parameter p defined by, 

p = i - (\+ +^_ ) 
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Here \ + and k_ are the eigenvalues of the matrix, 


U 


xy 




a 



For cellular two-dimensional flow p = 0 and for three dimensional isotropic 
flow p = 1 . 

The horizontal spectrum is generated by the following formula; 


a x N = NS , N = 1, 2, * • • 

a y M = M8 , M = 1, 2, • • • 

Here a x N anda y M are the x and y components of a. Truncation is achieved 
simply by discarding any a for |a| <a Q or |a| >a x . In the present experi- 
ment a Q = 1//2, a x = /2andS = 4 . The number of a_’s generated this way 

is 76 and the number of fluctuating-self-interactions treated is 6360 per vertical 
mode. 

In the present calculation the Rayleigh number is 3000, and 3 vertical wave 
numbers are used to describe the flow. This modest value of R, together with 
the small number of vertical wave numbers is dictated by machine time con- 
siderations. It may be objected at this point that R is too small to assess the 
accuracy of a statistical theory, since the flow is likely to be non-turbulent. Our 
point of view here is that even at these relatively small values of R a meaning- 
ful assessment of a statistical theory may be obtained by studying the relaxation 
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of the system from random initial conditions. A proper statistical theory should 
accurately treat the build up of correlations from an initial value of zero. We 
shall therefore insist that an acceptable statistical procedure should accurately 
treat the transient behavior of the system. 

In Fig. 11 and 12 we present results of the present numerical experiment 
forN u (t)and 0 (a, t). These results are an ensemble averaged over 10 

experiments. The initial value for the isotropy parameter is p - 0.88; at t = 

2.0, p - 0.80. The results are compared in the figures to results for the sta- 
tistical theories . For the latter we have not used the single band approximation 
to treat the horizontal wave number spectrum, but instead have used a more ac- 
curate treatment involving three a’ s . The procedure treats the a- integrations 
in (10a) and (10b) by the mid-point rule, and the a- integration in (2) by a Lagrange 
interpolation formula. The procedure is described in more detail in appendix C. 

The results for the direct interaction approximation for both N( t ) and 
4 n (a, t) appear to be satisfactory. The slight overestimation of N u at large t , 
together with a slight under estimation of \p 22 (a, t) may be traceable to the fact 
that the descretizing procedure for a omits an adequate representation of the 
linear harmonics (those terms for which a. = 2a.). These terms may play an 
important role in the energy transfer between 4> n and 0 22 , because these modes 
interact directly only through the mean field term. 

The quasi-linear approximation forN u (t) is also fairly good, especially 
considering the simplicity of this procedure. Its results for i p (a) for short 
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times is satisfactory but for long times it gives a0 n which sharpens into a 
delta function. Also, it predicts a </» (a) which tends to zero very rapidly as 

time increases (at t = 1.0, \p 22 ~ 10“ 6 , while \p — 10 -3 for either the numeri- 
cal experiment or the direct interaction). 

The quasi- normal approximation does not appear to behave satisfactorily; 
t/»n (a 3 , t) becomes negative for t - 0.3. 

The direct interaction results are sufficiently close to those of the numeri- 
cal experiment to make it of interest to inquire as to what significance should be 
attached to their difference. Some assessment of the numerical accuracy of the 
direct interaction integrations may be made by varying the method of integra- 
tion in the following two respects: (1) decrease the time step At, (2) use a dif- 
ferent integration method for the a-spectrum. With regard to the first point it 
appears that the time step used in the calculation is sufficiently small to guarantee 
good accuracy for the N u (t) curve: decreasing from 0.05 to 0.025 produces 
only a 0.01% difference in N u (t ) . Two tests were made to assess the accuracy 
of the a -integrations; first, the a -integration of Eq. (2) was modified from the 
Lagrange interpolation method to the mid point rule. This change produced a 
maximum difference inN u (t) of 0.26%. Next, the treatment of the fluctuating terms 
was changed by using a p (see Appendix C) averaged over Aa instead of an un- 
averaged one. The combined change produced by both modifications is at most 0.4%. 

The point to test in the numerical experiment is whether as the number of 
o-modes is systematically increased a significant change in any of the ensemble 
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quantities results. In order to test this point, a calculation was made using 
S = 1/5, which gives 124 a's distributed on (a Q < a < 2a 0 ). Some values of 
N u are indicated in Fig. 11 at the time for which the deviation of the 124 mode 
from the 76 mode experiment is maximum. The difference is at most 1.6%. 
However, it is not certain that this deviation is statistically significant since the 
standard deviation of N u in the 76 mode experiment is 2.2%. 

From the preceding paragraphs it appears that the maximum error for the 
direct interaction at R = 3x 10 3 is at most ^3%. The maximum occurs at 
t = 0.5, and the approximation appears better at larger times. 

VI. DISCUSSION OF RESULTS AND CONCLUDING COMMENTS 

A comparison of the DI results with the numerical experiments described 
in the last section suggests that the DI procedure is in good quantitative agree- 
ment with the exact solution to the thermal convection problem at large Prandtl 
number. The accuracy of the method on this point appears good for the transient 
behavior, as shown by Fig. 11 and 12, as well as for the steady state value. At 
R = 10 4 , the value of the Nusselt number predicted by the DI method is 5.03, as 
compared to 5.45 for the quasi-linear method, and 5.20 for the roll solution with 
a 0 = 1 // 2 . 

Although R = 10 4 is not considered a large value of Rayleigh number, it 
should be remembered that for the quasi-linear method it is sufficiently large 
for free boundaries for the dependence of N u on R to reach its asymptotic form. 
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This is also true for the two dimensional calculations of Fromm at a Prandtl 
number of unity. 

However R = 10 4 is not sufficiently large to draw conclusions about the asymp- 
totic shape of either the mean temperature field, or the temperature covariance 
field. The computed temperature profiles for R ^10 4 are properly bounded in the 
sense described in Section IV b, but violations of this bound if they occur at all 
would probably not appear until R = 10 s . 

The question of whether or not the numerical solutions of the Direct Inter- 
action equations are indeed static solutions is not entirely satisfactorily settled 
by the present calculation. It would be difficult to decide this issue entirely 
on the basis of numerical studies in any case, and the fact that the system 
evolves very slowly makes the problem that much more difficult. What is needed 
here is a stability analysis of the steady state solution at large Rayleigh numbers. 
Such a calculation seems rather formidable. Despite these uncertainties, we 
believe the numerical results suggest a non- static solution for the convection. 

The non-static character appears just above critical Rayleigh numbers, and 
becomes more pronounced as the Rayleigh number is increased. The occurrence 
of the time dependence apparently does not commence suddenly at a certain 
critical Rayleigh number as for example would occur in shear flow problems, 
if the Reynolds number is increased beyond a certain critical value. 

In Appendix B, it is pointed out that the DI approximation could give 
(static) results identical to the quasi-linear procedure provided R 1.18 x 10 4 . 
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The '/'-field would in these circumstances be static and non-statistical. Above 
R = 1.18 x 10 4 such a solution cannot be obtained, as is shown in Appendix 
B. The numerical results of the preceding section are again not conclusive, 
but they indicate that the final asymptotic state is not the quasi-linear type 
solution at any R. This conclusion is partly suggested by the shape of the 
correlation functions as given in Figs. 4, and 9. A comparison of 4 with 
9, suggests that as R increases, the solutions become less static in char- 
acter. The conclusion is reinforced by the fact that the numerical solution 
for R = 1.25 x 10 4 , which is above the critical value R = 1.12 * 10 4 beyond 
which the static quasi-linear solutions cannot exist show the same qualitative 
behavior as solutions below R = 1.12 * 10 4 . We therefore tentatively conclude 
that the quasi-linear-type solutions are not stable. 

The preceding paragraph tacitly assumes that the static solutions are the 
quasi-linear solutions, an assumption which is not necessarily true. In this con- 
nection, it may be pointed out that all efforts to find static solutions by iteration 
techniques always convergenced to the quasi-linear results. The attempt 
to find static solutions mentioned here used a vertical modes, one value of a, and 
R = 4 x 10 3 . 

Our investigation suggests that the infinite Prandtl number thermal con- 
vection problem is a peculiar limit case with some features are probably not 
to be found in the finite Prandtl number system at high Rayleigh numbers: 
namely, the extremely slow approach of the system to the steady state, and 
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the fact that the steady state solutions are nearly static. These features are not 
expected at finite Prandtl numbers, since the eddy viscosity then plays an im- 
portant role, and presumably introduces major decorrelation effects. The near 
static nature of the solutions suggests that the infinite Prandtl numbers regime 
may be an inhospitable one for statistical approximations. By the same token, 
it may therefore subject a statistical theory to a severe test. 

The analysis of the quasi-normal procedure indicates that this procedure 
is not an acceptable one to treat the evolution of nonstationary turbulence to 
the steady state for R £2 x 10 3 . The procedure may be an acceptable one to 
treat stationary turbulence for larger R, but the present results strongly suggest 
that if this is so, the method is very unstable to small perturbations and could 
not be used to treat small departures from equilibrium. 
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Appendix A 

Equations of Motion at Infinite Prandtl 
Number and Their Fourier Transform 

The equations determining the velocity field, v and temperature field T 
written in a convenient non-dimensional form are ; 

V • v = 0 

(^ 5 ^-V 2 )v = -Vp-I(v>-V)v-kRT 
[ji - V! ) T = - v ' VT 

These are the Boussinesq approximation to the Navier-Stokes equations. In 
this non-dimensional form the only constants which appear in these equations 
are the Prandtl number a - v/k and the Rayleigh number R = gaD 3 ATAv . Here 
a is the thermometric coefficient of expansion, v is the kinematic viscosity, k 
is the coefficient of thermometric diffusivity, AT is the temperature excess of 
the lower plate above the top one, D is the distance between the plates, and g is 
the gravitational acceleration. The non-dimensional variables v, T, r, t are 
related to the dimensional ones v' , T' , r ' , and t ' by; 

_ _ k _ 

v ~ D v 


(Al) 

(A2) 

(A3) 
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T' 


T 

T“ 

l o 


7 ' = Dr 


t' 


D 

IT t 


In the limit a -co, Eq. (A2) simplifies to a linear relationship between the velocity 
field and the temperature field, and thus in this limit v may be eliminated from 
the system. The elimination of the velocity field is facilitated by taking the curl 
of the curl of (A2). This leads to, 


V 4 v = RV (^ ) - R k V 2 T 


(A4) 


Slip boundary conditions for v are assumed on the conducting plates. Thus, 
if w is the vertical component of v: 


w 


<9 2 w 
dz 2 


= 0 


for 


and 


r - (x, y, 0) 


r - (x, y, 1) 
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The boundary conditions on the temperature field assume the boundary plates 
are infinitely conducting compared to the fluid they bound 

T(x, y, 0; t) = 0 

T(x, y, 1; t) = - 1 

Next, the temperature field T is split into its horizontal averages plus de- 
viation from this average. 


T = T + 6 

The velocity field has zero horizontal average . An average over the horizontal 
is indicated by the over bars. On using the above definition of 6 , Eq. (A3) splits 
into two components; 



/3w - V • ( v (9 - v<9 ) 


(A3a) 




(A3b) 


where ; 


fi s 


dT 

dz 


Now let the Fourier transform coefficients, w n “, 9 *, and T n be defined by 


w( r , t) 


L 

cl, n 


w a (t ) e l71a ' r sin nrrz 


(A4a) 
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6(r, t) - 


L e: 

cl ,n 


( t ) e 


Ittclt 


sin n 77 z 


(A5b) 


and 


T(z, t) 



T n (t)sinn7rz 


z 


(A5c) 


Here w( r , t ) is the vertical component of the velocity field. The horizontal 
component of v may be obtained in terms of 0 by using (A4) and (A5c) . The 
reality of v(r, t) andT(r, t) implies that 




and 


( T “)* = T" a . 

Here the asterisk denotes complex conjugation. The equations for w n a , 6 > n a , and 
T n may now be obtained by introducing (A5a), (A5b), and (A5c) into (A3a), and 
(A3b). The results are Eqs. (1), and (2) of the text. 
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Appendix B 


Connection Between Quasilinear and Direct Interaction Approximations 


a. Considerations Based on a Simple Model 

This section describes an analytic approach to the direct interaction equa- 
tions for a simple model system obtained by deleting all but the first two vertical 
Fourier modes and all but a single horizontal wave number a. This highly 
truncated system is investigated in the hope of gaining some insight into the 
structure of the full direct interaction equations. A direct analytic attack on 
either the quasi-normal or direct interaction equations is too difficult. The 

structure of the quasi-linear equations is sufficiently simple to yield to analytic 
4 

methods. The result of the present analysis is that for the simple two-mode 
model the static DI results for 0 nra are identical to the quasi-linear results. 
However, the two methods predict different results for g nm . 

Consider then the simple system for which the only vertical modes are 
n = 1, 2, so that the system (5), (10a), and (10b) consist of the eight functions 
0 n (t | t'), 0 12 (t | t'), 0 22 (t | t'), G n (t | t'), G 12 (t | t'), G 2 j ( t | t '), and 
G 22 (t | t ' ). The four off diagonal terms 0 12 , 0 21 , G 12> and G 21 are put equal to 
zero, since the statistical fields 0(7 | 7') may be assumed to be symmetric 
about the mid point, z = 1/2. The equations are, 


+v r) 0 ( t I t ' ) - Bn 0j ( t | t ' ) + M : 


j ’ 

•'n 


(t' I s) 0“ (t I s)0“(t I s) 


+ M" 


[ ds g“ < ( 
•'o 


f 

•'o 


(t I s) I s) i// a ( t I s) - NT ds g 2 a (t I s) \jJ x ( t ' I (t ! s) (Bl) 
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0 2 a (t I t') 


B 2 “ 2 ^2 (t I t') 


+ 2M 2 f ds g“ (t ' | s) 0“ ' (t | s) 0 2 a " (t | s) 


- 2M 2 I ds gj a ' (t | s) 0“ (t ' | s) 0“" (t | s) (B2) 
Jo 


(ar + v l °) gl a 


(t I t') = Bj a j g“ (t I t') 


- 2M 2 £ ( s | t')ds(g“' (t | s)0“" (t | s) 


+ g, a ' (t I s)0 2 a " (t | s>) (B3) 


{dt + ^2 a )g 2 a 


(t | t') = B 22 g“(t | t') 


2M 2 ds g“ (s | t')g 1 a '(t|s)0 1 a ”(t| s) (B4) 
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where, M = M U2 . In writing (Bl) - (B4), the condition M 121 +M 2U + M 212 = 0 
has been used. The a index is suppressed and only one vertical index is used to 
specify the G and 0 fields. The equations for 0( 1 1 t) are not recorded here; 
they are identical with (Bl) and (B2) with t = t' except that the term ( dp/d t ) 
is replaced by (1/2) (dp/ dt). 

It is now shown that in the limit t - °° a particular solution to (Bl) - 
(B4) consists of a static quasi-linear 0-field, for which 0 2 = 0, and <p x is specified 
by 


V 


a 

1 


B 


a 

1 1 


or 


01 = £ U ~ R c/ R ) - K = 657.7 . (B5) 

Examination of the steady state form of (Bl) - (B4) indicates that to prove 
this assertion it is necessary that 

02 ~ 0 (B6) 


and 



ds g 2 (s) 


0 . 


(B7) 
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Equation (B7) is a consequence of (B6) as may be seen by examining (B7) and 
using the assumption that as t - < p becomes a static field. 

To prove (B7) write the steady state form of (B3) and (B4) , and make a 
change of time scale, r = c 2 (t-t'), c 2 = 2M 2 </>! ; 


d 

d7 «i ( T > 



ds gj (r - s) g 2 (s) 


( 


d 

dF + v- 


) 


e 2 ( r ) 



dr g 2 (r - s) gj (s) 


where, 


* ~ (^2 ~ B 2 2 )/ c2 ’ 

These equations may be solved by using Laplace transform theory. Thus, 
introducing 


81 (s) 


g 2 ( s ) 


f.. 

•'O 


(t) e~ st dt 


/•CO 

Jo 


g 2 ( t ) e" st dt 


(B8) 


(B9) 


38 



algebraic equations for g t (s)and g 2 (s) maybe obtained, whose solutions are 


g x (s) 


2(s +m) 

|/s 2 ( s + /j,) 2 + 4s( s + fj.) + s( s + fj.) 


(BIO) 


g 2 ( s > " sg,(s)/(s + M) . (Bll) 

From (BIO) it follows that in the limit s - 0, g, ( s ) - Vs t and hence, from (B9) , 

s — 0 

(B7) is valid provided the limit process and Laplace transformation are 
interchangeable. 

Equations (BIO) and (Bll) may be inverted to give gj ( t ), and g 2 ( t ). It 
is of interest to have their asymptotic forms: 


g x (t) ^y£Ar " 1/2 (B12) 

7 — *00 

g 2 ( r ) * “ (l/2 it- 7 t)t~ 3/ 2 . (B13) 

7 * — *00 

This analysis assumes the stationary state system is static. As pointed out 
before, this assumption is mathematically consistent, but the static solutions 
may be unstable to small perturbations, which will induce the system to revert 
to a unstable state. 

The question now arises as to whether these simple results are genera- 
lizable to the complete multi-mode system. They may be. To see why, 
consider the DI equations assuming that the 0 -field is that found from 
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quasi-linear theory, < p h . It is known that this field has only a single a 
(say a 0 ), and that it produces (through Eq. (2)) a mean temperature distribution 
upon which any temperature fluctuation field other than *p L will decay. We shall 
now argue that in the limit t -*<» the assumption 4 ) ~ , 'P L is self consistent; that is 
to say, if <// L is inserted into the right hand side of the DI equations no new terms 
other than 0 L are generated at t -®. 

To establish this point, we need certain properties of the Green’s functions for 
\p - \p L , The particular properties of G n ^° (t | t') necessary are generalizations 
of the integrals of Gj andG 2 established for the simple model. These are, 


lim 

t-00 


lim 

t-^oo 



( t ' ) dt ' - +00 , 


(t') dt' - -0 , 


for (i, j) 


f°r (i, j) 


odd , (B14) 


odd . (B15) 


The even-odd components of g ; “ are zero, since i p L has no even-odd components. 
For a / a Qt the integrals are finite. 

The conditions under which (B14) and (B15) are valid are discussed in sec- 
tion b of the present appendix. The results derived there indicate that these equa- 
tions may be valid for R < 1.1 x 10 4 ; but cannot be correct for R > 1.1 x 10 4 . 

The consistency of the assertion that the DI solutions are 0 L is now examined. 
To do this, it will suffice to examine (Bl) - (B4) provided we recall that 
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these equations should have appropriate a integrals as given in (10a) and (10b). 
Note that in (Bl, B4)the fluctuating self interactions make two types of con- 
tributions to i p a ; the input term, containing the factor g a , and the drain term, 
containing the factor 0“ (t' | s). For example, in (Bl) the second term on the 
right hand side is an input term, and the last two terms are classified as drain 
terms. The drain terms may be positive (as is the third term on the right hand 
side of (Bl)) but in the main, they are negative. Similarly, the input terms are 
mainly positive (they are strictly positive if the non-diagonal elements of 0 and 

a 

g are suppressed). Consider first the odd-odd 0 equations. The input terms 
all vanish, since these contain even-even indexed 0 fields, and 0 L contains no 
such components. The drain term vanishes because of (B15) . Next consider 

CL 

the even-even 0 0 equation. The input term vanishes because of (B15) , and the 
drain term vanishes because there is no even-even component in 0 a ° . 

For a / a Q , the odd-odd components receive zero input, and a finite drain. 
Hence, they will vanish in the limit t -®>. The even-even components do receive 
an input, however, the drain term (the last term in (B2)) becomes indefinitely 
large because of (B14): hence the even-even components - 0 as t -<». The in- 
put terms cannot become indefinitely large because if a / a 0 the time integrals 
of G a are finite. 

The above discussion suggests that there may exist static solutions to the 
DI equations which in the limit t -> 00 which are identical with the quasi-linear so- 
lutions. This does not imply that if 0 L is introduced as initial data the solution 
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will remain equal to 0 L . This is because the DI equations involve the history of 
the system. The present discussion suggests that in this case, \p will develop 
transient values / \p L which will decay ast-®. 

Nothing has been said in the above discussion about the stability of the *p h 
field, or indeed of any static solution to the DI equations. Such information would 
be very valuable. 


b. Asymptotic Form of Green's Function for 

This section examines the form of the DI Green's function G. . (t - t ' ), for 
the quasi-linear intensities, </» L . The goal is to seek conditions under which (B14) 
and (B15) are valid. It is first shown that if G i . (t ) tends to a finite limit as 
t -oo, then 

G ii <‘> - . *i, > 0 

for 

( i , j ) = odd , (B 16 ) 

and 



(t)dt 


0 


for 


(i,j) - even . 

This is done by examining the Laplace transforms of G^ (t) : 


(B17) 


G . (s) 
ij v J 


■ r 
J 0 


e~ st dtG.j (t) 


(B18) 
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The behavior of G ij (t) for t -®> may then be inferred by means of the Abel 
theorem, 19 

lim - lim [g(s)S v+1 r(v+l)l (B19) 

t-+<» t v s- + o L J 

if v > -1.0, and 

lim t[G(t) -G(t - 1)1 - 0 , 
t -.00 L J 

and 

1 im e - *** G( t ) - 0 , 

t-o° 

for 

Real (p) > 0 

The Laplace transform of the direct interaction Green's function Eq. (10b) 
for i fi = \p , may be written in matrix form as, 

(s +^)G = I + CGCG ( B2 °) 

where ; 

n ~ - v S + B 

>n m nnm n m 

I =8 

nm nm 
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^nm ^ ' ^npm 

P 

here d p is the static quasi-linear temperature fluctuation amplitude. In deriving 
(B20) use of the condition \p nm = 9 n @ m for the quasi-linear system has been made. 
The mean field for the quasi-linear interaction, B nm is also made up of the 0 -field 
in accordance with (2). Recall now that 9 n has only odd components. Then by 
using the definition of M npq (see Eq. (2) ) it may be shown that C nm ¥ 0 only if 
n + m = odd. This permits (B18)to be split apart into even and odd components. 
Let P a be a projection matrix which deletes from any vector the even-Fourier 
components, and let P b be the compliment of P a , P a + P b = I. Then two equations 
for the even-even and odd -odd components of G may be found by pre and post- 


multiplying (B20) by first P a and then by P b : 


( sI . + O g . 

*a + ^ab ^ ba G a 

(B21) 

K'PbK 

*b + C ba ^a ^ab ^b 

(B22) 

where 



G 

a 

= P a GP a . 


i 

= P , 


a 



c ab 

= P. c Pb • 


c ba 

= PbCP, , etc. 
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The derivation of (B21) and (B22) uses the fact thatP a CP a = P b CP b = 0. Equations 
(B21) and (B22) may be collapsed by systematically deleting the zero interestitial 
rows and columns. The collapsed equations are identical in form to (B21 and (B22), 
with I a and I b replaced by I: 


(s+ Ma )G a = l + C ab G b C ba G a 

(B23) 

( S+ ^b) G b = 1 + C ba G a C ab G b 

(B24) 


The collapsed form of C's are 



The solution of (B23) and (B24) for G a is; 

G a = ( S+ ^a) _1 Q (B25) 

where 


Q 2 _ tq = - r 


and 


r = (»+M.)c b -> (s*M b )c,;> . 
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It is convenient to write the matrix elements of G a in the laminar eigenmode 
representation. The elements are 


lo.l i.) 


1 

2( st ^. k ) 


r (n ± fr? " 4 >i ) ( k . I n) (y, l i.) 


(B26) 


In Eq. (B26), /x a k is the k th eigenvalue of the odd-odd part of /x, | j a > is its as- 
sociated eigenvector, and< k a | denotes the eigenvector to the adjoint matrix, /x. 
The y ^s are eigenvalues of the matrix r, |y. > denotes the corresponding eigen- 
vector, and < y. \ the eigenvector of the adjoint matrix r . 

Equation (B24) maybe cast in alternate forms by using the properties of P; 


< k . lo.l j.> 1 2 1 lEl7 i)< 7 i |E " M i.> 


(B27) 


= ( 7 i * ' 4 ?i) < k a |El 7 i) ( 7 i 1 j a) 


(B28) 


where ; 


E = c b V(s + M b )c-> 

The sign of the radicals in (B26), and (B28) is to be chosen so that 

lim (< k l G | j >) “* 8.. , 

_ n \ a 1 a 1 J a / kj 

s u 
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which corresponds to the condition initial condition on G i . (t), G^ (0) = S. . 
Equations (B27) and (B28) may be established from (B26) by using the definition 
off and the relations, 


. . . ( s+ ^ a k ) . | . • 

= —7; ( k a l £ l n) 

< y i 1 j a) = 7TirT< 7 i |E_11 j a> 

We are now in a position to examine the limit of G a (s)ass-O. First, note 
from (B26) and (B27) that if y. , (k a |E| y.) , (y i 1 E _ 1 1 j a ) , (k a | y. ), and 
(y. I j a ) remain finite as s - 0, so does (k a | G a | j a ) for either /j.J or M a k ^ 0. 

We assume temporarily that the above is so and shall return later to a (numerical) 
verification. Thus, under the above conditions all the laminar eigenmode matrix 
elements of G will remain finite as s-0 with the exception of W = = 0 ). 

Furthermore, from (B25), it is expected that at least one of the 77 s * sa y J\ be- 
comes linear in s as s - 0, since the first factor in T, (f± a + s) has one zero 
eigenvalue as s - 0, namely that corresponding to /x k a = 0. We assume tem- 
porarily that y x "-as as s - 0, where a > 0, and suppose that y x (s ) remains 
negative and real for sufficiently small s . Then according to (B25) and (B26), 

W l G . I M.'> — lElr.Xn |E' 1 U.'> (B29) 
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Here (B28) is used only to exclude the possibility that G -* l/s as s- a>. 


Equations (B29) and (B19) may now be used to infer the behavior of G( t ) as 


t -oo: 


1 1 


W iG-l M, 1 } - Tv lElr.Xr, IE-U.') 


(B30) 


Or, in terms of the Fourier representation of G , 


G a / t N . __ A 

nm t-oo y^t nm 


(B31) 


where 


Equations (B30) and (B31) are valid provided the conditions of the Abel theorem 
are met; i.e., provided the Laplace transform integral of G n ( [ ^ ) ( t ) exists in the 
right half t -plain, exclusive of the imaginary axis. Examination of (B26) shows 
that a necessary condition for this to be so is that the quantity y i (s) {y i (s) - 4} 
does not become negative for any s > 0 and any i . If the above conditions are not 
satisfied, Eq. (B26) implies that G( s ) has a branch cut on the real s-axis which 
must intrude into the right half plane. The branch cut will make a contribution 
to G( t ) which will behave at large t like /t e at . In order to obtain sufficient 
conditions for (B31) to be valid, it must be shown that none of the possible sin- 
gularities of G(s) in the right hand plane contribute to the inverse Laplace 
transform. 
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Equation (B2) for the even-even components of may now be verified by 
inserting (B15) into (B9) and taking the limit as s -* 0. From this it follows that 

G ij ( s ) - /s B i3 , 

for 


( i , j ) - even . 

(B2) follows immediately by using (B3) in the limit s - 0. 

To complete the discussion of (Bl) and (B2) we must clear up several as- 
sumptions made in the discussion relating to the matrices A, E, and T. In order 
for (Bl) and (B2) to be valid, > 0, E must be non-singular, and the eigen- 
value spectrum of T(s) must be such as to make | ty? (s) - 4y. (s) real for 
s > 0. These points may be explored numerically in lieu of an analytic procedure. 

The numerical results on this point indicate that for R < 1.18 x 10 4 the as- 
sumptions enumerated in the above paragraph are valid. The computed y i spectra 
consists of negative real numbers, for sufficiently small s > 0, and explicit cal- 
culation shows that A. . > 0. ForR > 1.18 x 10 4 the behavior of y 1 (s) is not con- 
sistent with the above condition on the y^s. For R > 1.18 x 10 4 , y t is real and 
positive at small values of s and becomes real and negative as s - + °°. The 
behavior of y l (s ) for s-Oasa function of R is depicted in Fig. (Bl). 

The numerical calculations from which these conclusions were drawn consist 
of using an eight mode Fourier representation of the quasi-linear amplitude sys- 
tem, 6 n . The relevant matrices C ab , C ba , P, etc. were then truncated 4x4 
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matrices. Several runs using a four mode Fourier representation of 6 n indicated 
that the eight mode representation is adequate. 
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Appendix C 
Numerical Procedure 

This section describes the procedure for numerically integrating the direct 
interaction Eqs. (5), (10a), and (10b), starting from arbitrary initial conditions 
for the covariance matrix 0 nm (0 | 0). Integration of the equations for the quasi- 
linear and quasi-normal approximations may be made by simply setting to zero 
certain of the DI terms (as in the quasi-linear approximation) or replacing the 
certain terms by simplified expressions (as for example, in the quasi-normal 
method Eq. (11) is used to approximate the non-simultaneous correlation coef- 
ficients in terms of simultaneous ones). 

Equations (5), (10a), and (10b) are a set of integro-differential equations for 
^nm ( t I t ' ). They are discrete in the vertical index (n, m) and continuous in 
the indices (or arguments) a, and (t, t') . In treating these equations numerically, 
a suitable truncation procedure for the wave numbers (n, a), and a discretizing 
procedure for a, t, and t' must be prescribed. The procedure now described 
does this in a way so as to preserve the conservation properties of the flow. 

The flow is assumed to be isotropic in the horizontal direction, so that </< a 
and g“ depend only on the magnitude of a. 

To truncate the system in wave number (n, a) discard any ip nm or g nm func- 
tion if either n or m > M or if | a | > a 1 or | a | < a 0 , where a 0 < In practice, 
the finite machine storage capacity impose a rather severe limit on the number 
of individual dependent variables treated. In the present analysis it was necessary 
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to take M <5, to fit the program into the IBM 36065 fast storage. According to 
previous numerical results for the quasi -linear method, M = 5 is sufficient to 
accurately describe the mean temperature field up to a Rayleigh number of 
-R - 10 4 . 

To procede further, the basic Eqs. (5), (10a), and (10b) are reduced by using 
the assumption of horizontal isotropy. These are of the form 



B a i p a 


MJ 


da' da" 


a 


/?(a, a' , a" ) A(a, a' , a" ) g a \p a ' \p° 




p(a, a' , a" )C(a, a' , a" ) i// a g a ' >/' a " 

(Cl) 


+ v 2 ) g 2 - B a g a + J dsj J* - a ~ — p(a, a' , a")C(a, a' , a") g a < p a ‘ g a " 

(C2) 




(27raaa) J p a (a(n -p) p 


-'P 


a 

P . n+p 


) 


where ; 


<C3) 


p( a, a', a") - 4aa' a'y a. + a* + a" ) (a + a ' - a") (a + a" ~ a) (a' + a" - a) 

Here, the vertical indecese (n, m) as well as the time arguments are suppressed. 
The latter are (t ' | s), (t | s), and (t | s), in the order prescribed by (Cl), and 
(C2). The factor p(a, a , a") is the density of triangles whose sides are of 


52 


length a, a' , and a". It must be taken to be zero if any of the factors in its 
denominator is zero. 

a. Single Band Approximation to Horizontal Spectrum 
The horizontal wave numbers spectrum may be most easily approximated 
by assuming for a simple square wave form; 



0 , otherwise . (C4) 

The procedure used in this paper is as follows: 

Multiply (Cl) and (C2) by (2 na da) and integrate over a, using (C4) to ap- 
proximate the (//-spectrum in the convolution terms. 

Approximate g a (t | t' ) by a constant-in-a-value, g(t | t' ). There results 
equations for 0 and g which differ from (Cl) and (C2) by first deleting the a 
integrations (with the p-factor) and then using averaged coefficients in place of 
v a , B a , A andC. The average coefficients are, 


f a l 

= (vr(cLi 2 - cx 0 2 1 (2naaa)v n a = n 2 + -^(a 1 2 + a ( 2 ) 

Ja o 

m 2 ln~m 1 


J n - (R/V 4 )" 


In 


2 . 2 y 
n 4 + a i 

n 2 +a Q 2 


|(a 2 -a 0 2 ) - n 2 y/[(n 2 + a 0 2 ) (n 2 + a 2 )] 


A ~ SIS da da' da" p(a, a' , a") A(a, a' , a" )/jff p(a, a' , a" ) dada' da" 
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The mean field Eq. (C3) is similarly averaged with the result that J p replaces 
J p a . It may be verified directly that the above averaging procedure preserves 
the conservation property given by (3) . 

b. Approximation of a-lntegratio n Using Several 04 . 

In this case it is necessary to work with 0(a) and g(a) where, 

0(a) = 2rra0(a) 

and 

g(a) = 277a g(a) 

These quantities satisfy equations like (Cl) and (C2) except that p is replaced by 
p = p/(477 2 aa‘ a"). The initial condition on g is g(a, 0) = 2na. Next, (Cl) and 
(C2) are evaluated at a k = a Q + (2k - 1)A/N and the first integral in (Cl) is ap- 
proximated by, 


N 

A 2 2^ pK- V a k )A(a i , a., aj g(a.) 0(a.) 0<aJ 
j.k=l 

Here a Q and 2a 0 are lower and upper limits beyond which the integrand is con- 
sidered to make a negligible contribution to the integral. Similar approximations 
are made in the other integrals of (Cl) and (C2). 

The integration in (C3) is simple enough to make a more accurate integration 
formula practical. In this case we use the points 0^), 4>( a 2 )> anc * to ma k e 

a quadratic interpolation of 0(a). The effect of this is to change J n (a. ) to J n (i) 
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where. 


Jn(i> 


f 2a 0 ( a ~ a j) ( a ~ a k) 

Ja 0 ( a i ' a j )( a i- a k) Jn(a) 


These integrals are of standard form. 


( j . k f i) . 


c. Time Integration . 

The time integration of (Cl) and (C2) is straightforward and has been de- 
scribed elsewhere in more detail by Kraichnan. 19 The domain of the integration 
is t ^ t ' . For t = t',G(t|t')=l, and i/<(t | t ) is given by (Cl) with the term 
(d^p/dt ) replaces by (1/2) (di/</dt). Equations (Cl) involves </» n a m (t | t'), t' >t 
in the convolution integral in the first term on the right hand side. Such terms 
are evaluated by using the fact that i/> nm (t|t') = < p mn (t' 1 1 ), which follows from 
Eq. (4). 

Next, a mesh of points (t . , t . j is placed in the ( t , t ' ) plane, and, using 
them, the time integrals are approximated by the trapezoidal rule. The mesh 
point spacing A is taken to be uniform and the same in both time directions. To 
obtain | t . ) from <//(t ,_ 1 | t.j, formally integrate (Cl) over the interval 

( 1 i- i ’ t i ) » replacing the right hand side by its average value on ( t 1 , t . ) . The 


result is, 


( _ 

*Mt,) = e" v " A I ‘j ) * [ r (‘i) + F(ti-,)] ■ 


where F(t . ] represents the right hand side of (Cl). This is an implicit equation 
for 0 (t . | t since F(t . ) involves i/»(t . | t . ] itself. It must be solved by 
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iteration. The procedure used here is to approximate (t £ | t .) by the second 
order iteration using F(t . ) = F(t i _ 1 ) to initialize the right hand side. This 
guarantees accuracy of order A 2 . Aside from the exponential factors this 
scheme is the standard predictor - corrector procedure, in which the predictor 
is first order and the corrector is of second order. The presence of the ex- 
ponential factors stabilizes the behavior of the large wave number modes. 

The initial values of i p* m (0 | 0) are arbitrary. For most of the calculations 
they are chosen to be those values obtained from exact solutions to the hexagon 
amplitude equations, in which only a single horizontal wave number is included. 
That is to say. 


C. <° I 0) = ».«. <C5, 

This choice is a convenient one in that the system does not develop violent 
transient oscillations from these initial data. This permits a relatively large 
mesh size to be used. 

It is known that the hexagon solutions to the convection problem <9(7, t ), 
are asymmetric about the mid point z = 1/2, with the direction of asymetry ar- 
bitrary. This asymmetry is eliminated by averaging (C5) over the direction of 
asymmetry, which is equivalent to setting to naught the even -odd components of 
(C 5) . 
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FIGURE CAPTIONS 


Fig. 1 Nusselt number as a function of Rayleigh number R for: (a) the direct 
interaction, (b) quasi-linear, (c) quasi-normal methods. The calcu- 
lations contain 2 vertical and one horizontal degrees of freedom. 

Fig. 2 R n (8 | t)and R 22 (8 i t)for the direct interaction procedure for R = 10 3 . 

Fig. 3 Green's functions for direct interaction method for R = 10 3 . Param- 

eters of the calculation are the same as in Fig. 2 

Fig. 4 Direct interaction correlation coefficients R n (t| t-1.5) and 
R 22 (t | t - 1.5) for R = 4 x 10 3 . 

Fig. 5 Direct interaction correlation coefficients R 31 (t | t-1.5) and 
R 42 (t | t - 1.5) for same initial conditions as on Fig. 3. 

Fig. 6 Direct interaction Green's functions Gj 1 , G 13 , and G 31 for R = 4 x 10 3 . 

Fig. 7 Direct interaction Green's functions G 22 , and G 24 for R = 4 x 10 3 . 

Fig. 8 Evolution of Nusselt number for R = 10 4 according to Direct Interac- 
tion approximation. 

Fig. 9 Direct interaction correlation functions R n (t | 0.75 - t) and 
R 22 (t | 0.75- t) for R = 10 4 . 

Fig. 10 Evolution of < p n (t | t) according to the quasi -normal approximation 

at R = 2 x 10 3 . Run contains two vertical and one horizontal wave 
numbers. 

Fig. 11 Evolution of the Nusselt number for R = 3 x 10 3 , according to; (1) 

Numerical experiment (solid line), (2) The direct interaction (dashed 


57 



line), (3) Quasi-linear approximation points enclosed in triangles, 

(4) Quasi-normal approximation (dotted line) . Points enclosed in 
squares give N (t) for 124 a-mode experiment. 

Fig. 12 Evolution of ip lt ( a , t) for R = 3 x 10 4 . Solid curve represents the 
numerical experiment and the encircled points are direct interaction 
results. 

Fig. B1 (s)/s) for s = 10 -4 as a function ofR. 
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